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ABSTRACT 

In this paper, we compare and contrast the use of 
second-order response surface models and kriging models 
for approximating non-random, deterministic computer 
analyses. After reviewing the response surface method 
for constructing polynomial approximations, kriging is 
presented as an alternative approximation method for the 
design and analysis of computer experiments. Both 
methods are applied to the multidisciplinary design of 
an aerospike nozzle which consists of a computational 
fluid dynamics model and a finite-element model. Error 
analysis of the response surface and kriging models is 
performed along with a graphical comparison of the 
approximations, and four optimization problems are 
formulated and solved using both sets of approximation 
models. The second-order response surface models and 
kriging models — using a constant underlying global 
model and a Gaussian correlation function — yield 
comparable results. 

NOMENCLATURE 

(3 - constant underlying global portion of kriging model 
Pi, Py, Pii, - linear, interaction, and quadratic coefficients 
of polynomial equation in response surface 
DOE - design of experiments 
GLOW - gross lift off weight 
MSE - mean square error 
n s - number of sample points 
R - correlation matrix in kriging model 
R(x‘,x j ) - correlation function between points x 1 and x j 
RS - response surface 
6 2 - variance estimate 
0 k - correlation parameters in kriging model 
y - predicted response value at untried x 
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1 FRAME OF REFERENCE 

Current engineering analyses rely heavily on running 
expensive and complex computer codes. Despite the 
steady and continuing growth of computing power and 
speed, the complexity of these codes maintains pace 
with computing advances. Statistical techniques aie 
widely used in engineering design to construct 
approximations of these analysis codes; these 
approximations are then used in place of the actual 
analysis codes, offering the following benefits: 

• They yield insight into the relationship between 
output responses, y, and input design variables, x. 

• They provide fast analysis tools for optimization and 
design space exploration since the inexpensive-to-run 
approximations are used in lieu of the expensive-to- 
run computer analyses. 

• They facilitate the integration of discipline dependent 
analysis codes. 

A common method for building approximations of 
computer analyses is to apply design of experiments 
(DOE), response surface (RS) models, and regression 
analysis to build polynomial approximations of the 
computationally expensive analyses. For example, the 
Robust Concept Exploration Method 1 has been 
developed to facilitate quick evaluation of different 
design alternatives, identify important design drivers, 
and generate robust top-level design specifications using 
DOE, RS models, and the compromise Decision 
Support Problem; it has been successfully applied to 
the multiobjective design of a high speed civil 
transport, 2 a family of General Aviation aircraft, 3 a 
turbine lift engine, 4 and a flywheel. 5 In other work, the 
Variable Complexity Response Surface Modeling 
(VCRSM) method 6 uses analyses of varying fidelity to 
reduce the design space to the region of interest and 
build response surface models of increasing accuracy. 
The VCRSM method employs DOE and RS modeling 
techniques and has been successfully applied to the 
multidisciplinary wing design of a high speed civil 
transport, 7 to the analysis and design of composite 
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curved channel frames, 8 and to reduce numerical noise in 
staictural analyses. 9,10 A recent review of several 
applications of response surface models in engineering 
design is given in Ref. 11; for applications of 
approximations in structural design, see Ref. 12. 

Since computer experiments typically lack random 
error, a more appropriate and perhaps more “statistically 
sound” method for approximating deterministic 
computer experiments is through the use of kriging 
which is also referred to as the Design and Analysis of 
Computer Experiments (DACE) models. 13,14 The 
validity of the kriging model is not dependent on the 
existence of random error and may be better suited for 
applications involving computer experiments because it 
can either “honor the data,” providing an exact 
interpolation, or “smooth the data.” 15 

Booker 16 contrasts traditional DOE and RS 
modeling with DACE models. In the “classical” design 
and analysis of physical experiments, random variation 
is accounted for by spreading the sample points out in 
the design space and by taking multiple data points 
(replicates), see Figure 1. Sacks, et al. 14 state that the 
“classical” notions of blocking, replication, and 
randomization are irrelevant when it comes to 
deterministic computer experiments; thus, sample 
points should be chosen to fill the design space. 
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obtain output 


DOE/RS Modeling 
for Physical 
Experiments 


Account for Variability 



DACE/Kriging 
for Computer 
Experiments 

Space Filling 
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t-tests, F-statistics 
R 2 , R 2 adj ; Residual plots 
(see, e.g., Ref. 17) 


Maximum Likelihood 
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Cross-validation 
Integrated MSE 
(see, e.g.. Ref. 14) 


Figure 1. Comparison of DOE/RS Modeling 
and DACE/Kriging 16 


As noted in Figure 1, response surface modeling 
typically employs least squares regression to fit a 
polynomial model to the sampled data while kriging 
models are chosen to interpolate the data and are fit 
using maximum likelihood estimation. 13 Validation of 
RS models is based on: (a) testing statistical hypothesis 
(t-tests and F-statistics) derived from error estimates of 
the variability in the data, (b) plotting and checking the 
residuals, and (c) computing R 2 , the ratio of the model 
sum of squares to the total sum of squares, and R 2 adJ , 
which is R 2 adjusted for the number of parameters in the 


model. 17 Meanwhile, Sacks, et al. 14 and Welch, et al. 18 
both state that statistical testing is inappropriate when 
it comes to deterministic computer experiments which 
lack random error; therefore, cross-validation and 
integrated mean square error (MSE) are often employed 
to assess the accuracy of a kriging model. 

Kriging and DACE models have found limited 
application in engineering design perhaps because of the 
lack of readily available software to fit kriging models, 
the added complexity of fitting a kriging model, or the 
additional effort required to use a kriging model. To 
clarify this last point, prediction with a kriging model 
requires the inversion and multiplication of several 
matrices, and the kriging model does not exist as a 
“closed-form” polynomial equation; this is further 
clarified in Section 2.2. Meanwhile, RS model 
prediction requires computation of a simple polynomial 
equation once the model has been fit. The goal in this 
paper is to examine the added computational expense 
required to perform kriging and compare the predictive 
capability of kriging and RS models. 

In Section 2 an overview of the statistical and 
mathematical foundations of RS modeling and kriging 
is given. In Section 3 the multidisciplinary aerospike 
nozzle design example is introduced; it serves as a test 
problem to compare RS and kriging models for 
approximation. In Section 4 the RS and kriging 
models are constmcted and validated. In Section 5 four 
optimization problems are formulated and solved using 
the RS and kriging models, and Section 6 contains a 
discussion of ongoing work. 

2 STATISTICAL APPROXIMATIONS FOR 
COMPUTER EXPERIMENTS 

Building approximations of computer analyses typically 
involves: (a) choosing an experimental design to sample 
the computer analysis code, (b) choosing a model to 
represent the data, and (c) fitting the model to the 
observed data. There are a variety of options for each of 
these choices, and several of the advantages and 
disadvantages of each — with emphasis on response 
surface methodology, neural networks, inductive 
learning and kriging — are discussed in Ref. 1 1 . 

In this work, we are primarily concerned with the 
model choice and model fitting portion of building 
approximations. In particular, we focus on response 
surface models (Section 2.1) and kriging (Section 2.2). 

2.1 Overview of Response Surface Modeling 

2.1.1 Mathematics of Response Surface Modeling 
Response surface modeling techniques were originally 
developed to analyze the results of physical experiments 
and create empirically-based models of the observed 
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response values. Response surface modeling postulates 
a model of the form: 


yfx) = f(x) + e (1) 


where y(x) is the unknown function of interest, ffx) is a 
known polynomial function of x, and e is random error 
which is assumed to be normally distributed with mean 
zero and variance cr. The individual errors, £;, at each 
observation are also assumed to be independent and 
identically distributed. The polynomial function, f(x), 
used to approximate y(x) is typically a low order 
polynomial which is assumed to be either linear, Eqn. 
(2), or quadratic, Eqn. (3). 

y = P 0 +t(3 lXl (2) 

i=l 


' = Po + £ M. + £ Mi 2 + X E PijXjXj 


(3) 


The parameters, (3 0 , (3,, [3 H , and (3^, of the 
polynomials in Eqns. (2) and (3) are determined through 
least squares regression which minimizes the sum of the 
squares of the deviations of predicted values, y (x), from 
the actual values, y(x). The coefficients of Eqns. (2) 
and (3) used to fit the model can be found using least 
squares regression given by Eqn. (4): 

P=[X’X]' 1 X’y (4) 

where X is the design matrix of sample data points, X’ 
is its transpose, and y is a column vector containing the 
values of the response at each sample point. For more 
details on least squares regression or polynomial RS 
modeling see, e.g., Ref. 17. 

2.1.2 General RS Modeling Approach The general 
approach for building polynomial response surface 
models is shown in Figure 2. A three step process 
involving screening, model building, and model 
exercising is typically employed. 

As shown in Figure 2, the first step involves 
screening which may be employed if there are a large 
number of factors to reduce the design space to an 
appropriate region of interest. In the second step, the 
approximation models are built from sample data which 
is obtained from an appropriately chosen experimental 
design; if there are noise factors in the design, 
robustness models of the mean and variance of each 
response would also be created. If the models are 
sufficiently accurate, the model is exercised in the last 
stage of the process to search the design space and find 
improved or robust solutions; the reader is referred to, 
e.g., Ref. 17 for more information on response surface 
modeling. 


Model 

Building 



Model 

Exercising 


Given : 

Factors, Responses 



Solutions 

(improved or robust) 


Screening 


Figure 2. General RS Modeling Approach 11 


2.2 Overview of Kriging 

2.2. 1 Mathematics of Kriirinp Kriging postulates 
a combination of a global model plus departures: 

y(x) = f(x) + Zfx) (5) 

where yfx ) is the unknown function of interest, ffx) is a 
known (usually polynomial) function of x, and Zfx) is 
the realization of a stochastic process with mean zero, 
variance cr, and non-zero covariance. The ffx) term in 
Eqn. (5) is similar to the polynomial model in a 
response surface and provides a “global” model of the 
design space. In many cases ffx) is simply taken to be 
a constant term p (see, e.g., Refs. 13, 14, 18); we use a 
constant term for ffx) in the example in Section 4. 

While ffx) “globally” approximates the design 
space, Zfx) creates “localized” deviations so that the 
kriging model interpolates the n s sampled data points. 
The covariance matrix of Zfx) is given by Eqn. (6). 

CovfZfx^Zfx 1 )] = cr R([R(x‘,x J )] (6) 

In Eqn. (6), R is the correlation matrix, and Rfx'.x 1 ) is 
the correlation function between any two of the n, 
sampled data points x 1 and x j . R is a fn s x nj 
symmetric matrix with ones along the diagonal. The 
correlation function R(x‘,x j ) is specified by the user; 
Sacks, et al. 14 and Koehler and Owen 13 discuss several 
correlation functions which may be used. In this work, 
we employ a Gaussian correlation function of the form: 

R(x\x ] ) = exp[-r k :e k K “ x k|" ] ( 7 ) 
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where n dv is the number of design variables, 0 k are the 
unknown correlation parameters used to fit the model, 
and the x k ‘ and x k j are the k* 11 components of sample 
points x 1 and x j . In some cases using a single 
correlation parameter gives sufficiently good results 
(see, e.g., Refs. 19, 20, and 14). 

Predicted estimates, y(x), of the response y(x) at 
untried values of x are given by: 

y = P + r T (x)R- 1 (y-fp) (8) 

where y is the column vector of length n s which 
contains the sample values of the response, and f is a 
column vector of length a, which is filled with ones 
when f(x) is taken as a constant. In Eqn. (8), r T (x) is 
the correlation vector of length n s between an untried x 
and the sampled data points (x 1 , ..., x" 4 5 }: 

r T (x) = [R(x,x'), R(x,x 2 ), ..., R(x,x ns )] T . (9) 
In Eqn. (8), p is estimated using Eqn. (10). 

P = (f T R~ 1 f)” 1 f T R“ 1 y (10) 

The estimate of the variance, o 2 , between the 
underlying global model p and y, is estimated as: 

p (y-f(3)' r R -’(y-ip.) (]1) 

n s 

where f(x) is assumed to be the constant p. The 
maximum likelihood estimates (i.e., “best guesses”) for 
the 0 k in Eqn. (7) used to fit the model are found by 
maximizing Eqn. (12) over 0 k > 0 (see, e.g., Ref. 19). 

[n s ln(o 2 ) + In | R |] ^ 

Both d 2 and |R| are functions of 0 k . While any values 
for the 0 k create an interpolative model, the “best” 
kriging model is found by solving the k-dimensional 
unconstrained non-linear optimization problem given by 
Eqn. (12). 

2.2.2 Engineering Applications of Kriging DACE 
and kriging models have found limited use in 
engineering design applications since its introduction 
into the literature by Sacks, et al. 14 Giunta 21 performs a 
preliminary investigation into the use of DACE 
modeling for the multidisciplinary design optimization 
of a High Speed Civil Transport aircraft. He explores a 
5 and a 1 0 variable design problem, observing that the 
DACE and response surface modeling approaches yield 
similar results due to the quadratic trend of the 
responses. Booker, et al. 22 solve a 31 variable 
helicopter rotor structural design problem; Booker 16 
expands the problem to include 56 variables to examine 
the aeroelastic and dynamic response of the rotor. Osio 


and Amon 20 use a multistage DACE modeling strategy 
to design an embedded electronic package which has 5 
design variables. Finally, some researchers have 
employed DACE modeling strategies specifically for 
numerical optimization (see, e.g., Refs. 23 and 24). 

3 AEROSPIKE NOZZLE EXAMPLE 

The multidisciplinary design of an aerospike nozzle has 
been selected as the test problem for comparing the 
predictive capability of RS and kriging models. The 
linear aerospike rocket engine is the propulsion system 
proposed for the VentureStar 25 Reusable Launch Vehicle 
(RLV) which is illustrated in Figure 3. 



Figure 3. VentureStar RLV with Linear 
Aerospike Propulsion System 26 


The aerospike rocket engine consists of a rocket 
thruster, cowl, aerospike nozzle, and plug base regions 
as shown in Figure 4. The aerospike nozzle is a 
truncated spike or plug nozzle that adjusts to the 
ambient pressure and integrates well with launch 
vehicles. 27 The flow field structure changes dramatically 
from low altitude to high altitude on the spike surface 
and in the base flow region. 28-30 Additional flow is 
injected in the base region to create an aerodynamic 
spike 31 which gives the aerospike nozzle its name and 
increases the base pressure and contribution of the base 
region to the aerospike thaist. 



Figure 4. Aerospike Components and Flow 
Field Characteristics 26 
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The analysis of the aerospike nozzle involves two 
disciplines: aerodynamics and structures; there is an 
interaction between the structural displacements of the 
nozzle surface and the pressures caused by the varying 
aerodynamic effects. Thrust and nozzle wall pressure 
calculations are made using computational fluid 
dynamics analysis and are linked to a structural finite- 
element analysis model for determining nozzle weight 
and structural integrity. A mission average engine 
specific impulse and engine thrust/weight ratio are 
calculated and used to estimate vehicle gross-lift-off- 
weight (GLOW). The corresponding multidisciplinary 
domain decomposition is illustrated in Figure 5. Korte, 
et al. 26 provide additional details on the aerodynamic and 
structural analyses for the aerospike nozzle. 


GLOW contours 



Figure 5. Multidisciplinary Domain 
Decomposition 26 


For this study, we consider three design variables 
for the multidisciplinary design of the aerospike nozzle: 
(starting) thruster angle, (exit) height, and length as 
shown in Figure 6. 


nozzle surface profile based on the values of thruster 
angle, height, and length. 

Bounds for the design variables are set to produce 
viable nozzle profiles from the quadratic model based on 
all combinations of thruster angle, height, and length 
within the design space. Second-order response surface 
models and kriging models are developed for each 
response (thrust, weight, and GLOW) in the next 
section; optimization of the aerospike nozzle using the 
RS and kriging models for different objective functions 
is performed in Section 5. 

4 APPROXIMATIONS FOR THE 
AEROSPIKE NOZZLE PROBLEM 

The data used to fit the RS and kriging models is 
obtained from a 25 point random orthogonal array; 32 
The use of these orthogonal arrays is based, in part, on 
the work in Ref. 19 and the discussion in Ref. 33. The 
sample points are illustrated in Figure 7 and are scaled 
to fit the design space defined by the bounds on the 
thruster angle (a) , height (h), and length (1). 




The thruster angle is the entrance angle of the gas 
from the combustion chamber onto the nozzle surface; 
the height and length refer to the solid portion of the 
nozzle itself. A quadratic curve defines the aerospike 


In Section 4. 1 we discuss the RS models which ate 
fit to the data and in Section 4.2, the kriging models. 
Error analysis of the RS and kriging models is discussed 
in Section 4.3, and a graphical comparison of the 
approximations is given in Section 4.4. 

4.1 Response Surface Models 

Second-order RS models for weight, thrust, and GLOW 
are obtained using ordinary least squares regression 
techniques. 34 The corresponding RS models are given 
in the Eqns. (13)-(15). The equations have been scaled 
against the baseline design to protect the proprietary 
nature of some of the data. 

Weight = 0.810 - 0.116a + 0.121h (13) 

+ 0.1521 + 0.065a 2 - 0.025ah + 0.0013h 2 
- 0.0539al - 0.0131 111 + 0.03011 2 

Thrust = 0.9968 + 0.00031a + 0.0019h (14) 

+ 0.00601 - 0.00175a 2 + 0.00125ah - 0.001 lh 2 
+ 0.00125al - 0.00198hl - 0.001651 2 
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GLOW = 0.9930 - 0.0270a + 0.0065h (15) 

- 0.02651 + 0.0307a 2 - 0.0163ah + O.OlOOh 2 
- 0.0226al + 0.0151 111 + 0.01951 2 


The resulting R 2 , R 2 adj , and root MSE values for each 
RS model are given in Table 2. The root MSE is: 


root MSE = 


I^y.-y ,) 2 

i 


(16) 


where n is the number of sample points, y ; is the actual 
value of the response, and y ; is the predicted value. As 
evidenced by the high R 2 and R 2 adj values and low root 
MSE values, the second-order RS models appear to 
capture a large portion of the observed variance. 


Table 2. Model Diagnostics of RS Models 


Response 


Measure 

Weight 

Thrust 

GLOW 

R 2 

0.986 

0.998 

0.971 

R 2 aJ] 

0.977 

0.996 

0.953 

root MSE 

1.12% 

0.01% 

0.25% 


4.2 Kriging Models for the Aerospike 
Nozzle Problem 

Kriging models are built from the same 25 sample 
points used to fit the response surface models in Section 
4.1. We chose to model the data using a constant term 
for the global model and a Gaussian correlation 
function, Eqn. (7), for the local departures determined by 
the correlation matrix, R. 

Initial investigations revealed that a single 9 
parameter was insufficient to accurately model the data 
due to scaling of the design variables (a similar problem 
is discussed in Ref. 2 1 ). Therefore, an exhaustive grid 
search with a refinable step size was used to find the 
maximum likelihood estimates for the three 9 
parameters needed to obtain the “best” kriging model. 
The resulting maximum likelihood estimates for three 9 
parameters for the weight, thrust, and GLOW models 
are summarized in Table 3; these values are for the 
scaled sample points. 


4.3 F.rror Analysis of Response Surface and 
Kriging Models 

An additional 25 randomly selected validation points are 
used to verify the accuracy of the RS and kriging 
models. Error is defined as the difference between the 
actual response from the computer analysis and the 
predicted value from the RS or kriging model. The 
maximum absolute error, the average absolute error, and 
the root MSE — from Eqn. (16) where n (= 25) is the 
number of validation points — for the 25 randomly 
selected validation points are summarized in Table 4. 

Table 4. Error Analysis of Approximations 


Second-Order Response Surface Models 



W eight 

Thrust 

GLOW 

Max ABS(error) 

19.57% 

0.032% 

3.68% 

Avg ABS(error) 

2.44% 

0.012% 

0.53% 

root MSE 

4.54% 

0.015% 

0.90% 


Kriging Models (w/Constant and Gaussian Cor. Fen.) 



Weight 

Thrust 

GLOW 

Max ABS(error) 

17.23% 

0.048% 

3.43% 

Avg ABS(error) 

2.51% 

0.012% 

0.59% 

root MSE 

4.37% 

0.018% 

0.89% 


For the weight and GLOW responses, the kriging 
models have lower maximum absolute errors and lower 
root MSEs than the RS models; however, the average 
absolute error is slightly larger for the kriging models. 
As for thrust, the RS models are slightly better than the 
kriging models according to the values in the table; the 
maximum absolute error and root MSE are slightly less 
while the average absolute errors are essentially the 
same. It is not surprising that the RS model predicts 
thrust better; it has a very high R 2 value (0.998) and 
low root MSE (0.01%). It is reassuring to note, 
however, that the kriging model, despite using a 
constant term and a Gaussian correlation function, is 
only slightly less accurate than the corresponding RS 
model. In summary, it appears that both models predict 
well with the kriging models having a slight advantage 
in accuracy because of the lower root MSE values. 


Table 3. 9 Parameters for Kriging Models 


Response 

Weight 

Thrust 

GLOW 

0.548 

0.30 

3.362 

1.323 

0.50 

2.437 

2.718 

0.65 

0.537 


With these parameters for the Gaussian correlation 
function and the 25 sample data points, the kriging 
models are fully specified. A new point is predicted 
using these values in combination with Eqns. (B)-(10). 


4.4 Graphical Comparison of RS and 
Kriging Models 

In addition to the error analysis of Section 4.3, a 
graphical comparison of the RS and kriging models is 
performed to visualize differences in the two 
approximations. In Figures 8-11, three-dimensional 
contour plots of thrust, weight, and GLOW as a 
function of angle, length, and base height are given. In 
each figure, the same contour levels are used for the RS 
and kriging models so that the shapes of the contours 
can be compared directly. 
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Figure 8. Thrust Approximation Contours 


In Figure 8, the contours of thrust for the RS and 
kriging models are very similar. As evidenced by the 
high R 2 and low root MSE values, we expect the RS 
models to fit the data quite well. It is again reassuring 
to note that the kriging models resemble the RS models 
even through the underlying global model for the 
kriging models is just a constant term. 


Figure 9. Weight Approximation Contours 


The contours of the RS and kriging models in 
Figure 9 are also very similar, but we begin to see 
localized perturbations caused by the Gaussian 
correlation function in the kriging model for weight. 
The error analysis from the previous section indicated 
that the kriging model for weight is slightly more 
accurate than the RS model which might be attributed 
to these small non-linear localized variations. 


VicfeM 


Ot*r 

K' '.I.. . I 


Figure 10. GLOW Approximation Contours 


The general shape of the GLOW contours is the 
same in Figure 10; however, the size and shape of the 
different contours, particularly along the length axis, are 
quite different. The end view along the length axis in 
Figure 1 1 further highlights the differences between the 
two models. Notice also in Figure 1 1 that the kriging 
model predicts a minimum GLOW located within the 
design space centered around Height = -0.8, Angle = 0, 
along the axis defined by 0.2 < Length < 0.8; this 
minimum was verified through additional experiments. 


CLOU 


2ml Dnlf r 
HS Sfedtl 


■«gir 

Figure 11. GLOW Approximation 
Contours — End View 

The true test of the accuracy of the RS and kriging 
models comes when the approximations are used during 
optimization. This is performed in the next section. 
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5 OPTIMIZATION USING RESPONSE 
SURFACE AND KRIGING MODELS 

It is paramount that any approximations used in 
optimization are reasonably accurate, lest they lead the 
optimization algorithm into regions of bad designs. 
Trust Region approaches (see, e.g., Ref. 35) and the 
Model Management framework (see e.g., Refs. 36 and 
37) are being developed to ensure that optimization 
algorithms are not led astray by inaccurate 
approximations. In this work the focus has been on 
developing the approximation models, particularly the 
kriging models, and not on the optimization itself. 

We formulate and solve four different optimization 
problems to compare the accuracy of the RS and kriging 
models: (1) maximize thrust, (2) minimize weight, (3) 
minimize GLOW, and (4) maximize thrust/weight ratio. 
The first two objective functions represent traditional 
single objective, single discipline optimization 
problems. The second two objective functions are more 
characteristic of multidisciplinary optimization; 
minimizing GLOW or maximizing the thrust/weight 
ratio requires trade-offs between the aerodynamics and 
structures disciplines. For each objective function. 


constraint limits are placed on the remaining responses; 
for instance, constraints are placed on the maximum 
allowable weight and GLOW and the minimum 
allowable thrust/weight ratio when maximizing thrust. 
However, none of the constraints are active in any of 
the final results. The optimization results are 
summarized in Table 5. 

As shown in Table 5, each optimization problem is 
solved using: (a) the RS model approximations and (b) 
the kriging model approximations for thrust, weight, 
and GLOW. The optimization is performed using the 
Generalized Reduced Gradient (GRG) algorithm in 
OptdesX. 38 Three different starting points are used for 
each objective function (the lower, middle, and upper 
bounds of the design variables) to assess the average 
number of analysis and gradient calls to the 
approximations that is necessary to obtain the optimum 
design within the given design space. The same 
parameters (i.e., step size, constraint violation, etc.) are 
used within the GRG algorithm for each optimization. 
Design variable and response values have been scaled as 
a percentage of the baseline design to protect the 
proprietary nature of some of the data. 


Table 5. Optimization Results using Response Surface and Kriging Models 



Avg. # of 

Avg. # of 



Verified 



Analysis 

Calls 

Gradient 

Calls 

Optimum Design 

Predicted Optimum 

Optimum 1 

% Error* 


Maximize Thrust 


RS 

Models 

27 

4 

Angle 

Height 

Length 

™^Vrrrl 

Thrust 

Weight 

Thr/Wt 

GLOW 

1.0016 

0.9450 

1.0141 

0.9724 

1.0013 

0.9476 

1.0134 

0.9759 

in 




Angle 


Thrust 

1.0016 

1.0014 


Kriging 

62 

5 

Height 


Weight 

0.9385 

0.9155 

2.51% 

Models 



Length 


Thr/Wt 

1.0157 

1.0210 








GLOW 

0.9690 

0.9683 

0.08% 


Minimize Weight 


RS 

Models 

29 

3 

Angle 

Height 

Length 


Thrust 

Weight 

Thr/Wt 

GLOW 

0.9957 

0.7584 

1.0533 

0.9936 

0.9957 

0.7496 

1.0555 

0.9906 

-0.01% 

1.18% 

-0.21% 

0.30% 




Angle 


Thrust 

0.9965 

0.9956 

0.08% 

Kriging 

43 

4.67 

Height 

-0.873 

Weight 

0.7725 

0.7443 

3.79% 

Models 



Length 

- 1.000 

Thr/Wt 

1.0506 

1.0568 

-0.59% 






GLOW 

0.9824 

0.9914 

-0.90% 


Minimize GLOW 


RS 

Models 

30.67 

3.33 

Angle 

Height 

Length 


Thrust 

Weight 

Thr/Wt 

GLOW 

1.0013 

0.8969 

1.0251 

0.9660 

0.9957 

0.8617 

1.0286 

1.0146 

6.56% 

4.09% 

-0.34% 

-4.79% 






Thrust 

1.0009 

1.0006 

0.04% 

Kriging 

57.67 

6.33 

Hi SHmPM 

•SSk 

Weight 

0.9060 

0.8732 

3.75% 

Models 



sum 


Thr/Wt 

1.0228 

1.0302 

-0.72% 




■MB 


GLOW 

0.9675 

0.9680 

-0.05% 


Maximize Thrust/Weight Ratio 


RS 

Models 

27 

4 

Angle 

Height 

Length 

0.096 

-0.433 

1.000 

thrust 

Weight 

Thr/Wt 

GLOW 

1.0016 

0.9450 

1.0141 

0.9724 

0.9959 

0.9073 

1.0173 

1.0228 

0.57% 

4.16% 

-0.31% 

-4.93% 




Angle 

0.656 

Thrust 

1.0016 

1.0014 

0.02% 

Kriging 

62 

5 

Height 

-0.627 

Weight 

0.9385 

0.9063 

3.56% 

Models 



Length 

1.000 

Thr/Wt 

1.0157 

1.0231 

-0.73% 






GLOW 

0.9690 

0.9666 

0.25% 


^The predicted optimum value is obtained by using the values of angle, height, and length from the optimum design) in the actual analysis code. 
# A (+) error term indicates that the model is over-predicting; a (-) indicates under-predicting. 
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The following observations are made based on the 
data in Table 5. 

• Average number of analysis and gradient calls: In 

general, the optimization requires fewer analysis and 
gradient calls to the RS models than the kriging 
models. This can be attributed to the fact that the RS 
models are simple second-order polynomials whereas 
the kriging models are more complex, non-linear 
functions. 

• Convergence rates: Although not shown in the table, 
optimization using the RS models tends to converge 
more quickly than when using kriging models. This 
can be inferred from the number of gradient calls 
which is one to three calls fewer for the RS models 
than the kriging models. 

• Optimum designs: The optimum designs obtained 
from the RS and kriging models are essentially the 
same for each objective function. The largest 
discrepancy is the length when minimizing GLOW; 
RS models predict the optimum GLOW occurs at the 
upper bound on length (+1) while the kriging models 
yield 0.676. This difference is evident in Figures 10 
and 1 1 . 

• Predicted optima and prediction errors: To check the 
accuracy of the predicted optima, the optimum design 
values for angle, height, and length are used as inputs 
into the original analysis codes and the percentage 
difference between the actual and predicted values is 
computed. The prediction error is less than 5% for all 
cases and is 0.5% or less in three quarters of the 
results. 

In summary, the RS and kriging approximations 
yield comparable results with minimal difference in 
predictive capability. It is worth noting that the kriging 
models perform as well as the second-order RS models 
even though the global portion of the kriging model is 
only a constant. Ongoing work to further improve the 
accuracy of the kriging models is discussed next. 

6 CLOSING REMARKS 

This simple, yet realistic, engineering example of the 
design of an aerospike nozzle has been utilized to 
demonstrate the use of kriging models as an alternative 
approximation technique to second-order response 
surface models for multidisciplinary design 
optimization. There are several research issues to 
address for the application of kriging to other (and 
larger) engineering design problems. 

• Selecting a kriging model: In this example, we use a 
constant for the global portion of the kriging model. 
However, using a global, low-order polynomial 
model for f(x) in Eqn. (5) may further improve the 
accuracy of the kriging model. Giunta 21 performs a 


preliminary investigation of such an approach; his 
results indicate that minimal improvement in model 
accuracy is obtained. 

• Fitting a kriging model: For small problems with 
relatively few sample points, fitting a kriging model 
by optimizing Eqn. (12) is rather trivial. However, 
as the size of the problem increases and the number of 
sample points increases, the added effort needed to 
obtain the “best” kriging model may begin to 
outweigh the benefit of building the approximation. 

• Predicting with a kriging model: Unlike RS model 
prediction, prediction with a kriging model requires 
the inversion and multiplication of several matrices 
which grow with the number of sample points. 
Hence, for large problems prediction with a kriging 
model may become computationally expensive as 
well. The kriging software we are developing will 
facilitate fitting, building, and validating kriging 
models, increasing their attractiveness for engineering 
applications. 

• Validating a kriging model: Since kriging models 
interpolate the data, R 2 values and residual plots 
cannot be used to assess model accuracy. In this 
example we use an additional 25 random validation 
points to check model accuracy; however, other 
approaches exist. Yesilyurt and Patera 39 and Otto, et 
al. 40 have developed a Bayesian-validated surrogate 
approach which systematically uses additional 
validation points to make quantitative assessments of 
the quality of the approximation model and provide 
theoretical bounds on the largest discrepancy between 
the model and the actual computer analysis. An 
alternative method which does not require additional 
points is leave-one-out cross validation, 41 but it is 
uncertain how well this measure correlates with 
model accuracy. 

• Design of experiments for building kriging models: 
As discussed in Section 1, “space filling” 
experimental designs may be better suited for 
computer experiments. In this example, we use 
orthogonal arrays, but several other experimental 
designs exist. Koehler and Owen 13 discuss a wide 
variety of designs including Latin hypercubes, 
minimax/maximin designs and orthogonal arrays. 

Future work on the aerospike nozzle design 
problem includes adding more design variables and 
responses and investigating the impact of decomposing 
the problem into its disciplines by building separate 
approximations of each discipline and examining the 
effects of different multidisciplinary design formulations 
on the solution. 
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